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INTRODUCTION 


This report summarizes progress made on NASA Grant NAG-1-392 during 
the period February 19, 1986 through August 18, 1986. The goal of the 

project is to develop mathematical models and associated solution 
procedures which can be used to design heat pipe cooled structures for 
use on hypersonic vehicles. The models should also have the capability 

to predict off-design performance for a variety of operating conditions. 
It is expected that the resulting models can be used to predict startup 
behavior of liquid metal heat pipes to be used in reentry vehicles, 
hypersonic aircraft and space nuclear reactors. 

Previous reports in this series (February 1985, September 1985 and 
June 1986) have covered the development of governing differential 
equations which can be used to predict startup behavior of initially 
frozen liquid metal heat pipes. The June 1986 report covered the 
development of finite element based numerical equations for the outer 
shell and the combination capillary structure and working fluid. One and 
two dimension transient solutions were presented. 

The current report covers mathematical models and sample solutions 
for the vapor region. The heat pipe considered is rectangular, the 
working vapor is sodium and the initial temperature is sufficiently low 
so that the capillary structure is filled with solid sodium. 

ANALYSIS OF VAPOR FLOW IN HEAT PIPES 
Analysis of the hydrodynamics of vapor flow in heat pipes with 
metallic working fluids indicates considerable difficulties at low vapor 
pressure due to the extremely small vapor densities. Even for 
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relatively small heat transfer rates, the vapor velocity in the axial 
direction can be very large and accelerated towards sonic velocity by 
friction at the interface. Additional heat input causes choking at 
the condenser inlet. Such behavior makes it necessary to include vapor 
compressibility and viscous action in mathematical models. Vapor 
pressure drop due to friction can not be recovered completely in the 
condensation section. Thus, the temperature distributions along 
the length of the heat pipe are not isothermal and thermal resistance 
in the vapor region is significant. Therefore, studies of the 
distributions of temperature, pressure, and velocity in the vapor 
passage along the length of a heat pipe are essential for an evaluation 
of the maximum heat transfer rates and prediction of correct 
heat pipe performance. 

A summary of the literature on hydrodynamic processes in vapor flow 
of cylindrical heat pipes made by Tien [1] and Ivanovskii et. al . [2] 
indicates that no completely detailed investigation has been 
presented thus far. Moreover, experimental measurements of pressure 
and velocity for metallic working fluids have not been reported. 

In the present paper, an analysis of the steady, 
compressible, one-dimensional, laminar flow of sodium vapor is presented 
for the case of a flat plate-type heat pipe with asymmetrical 
boundary conditions. In addition, shear stress at the liquid-vapor 
interface, variations of vapor quality, and momentum and energy 
factors are considered. A similarity solution for a semiporous 
channel is used to provide the velocity profile at cross sections. 
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1. Similarity Solution for Serai porous Channel 
Governing Differential Equations 

A sketch of the geometry fora semiporous channel is given in Figure 
1. The following assumptions are made : that the fluid flow is 

incompressible and laminar, and that the properties of the fluid are 
constant. The width of the channel is assumed much greater than 
the height. Therefore, two-dimensional flow is considered. Fully 
developed flow is assumed in the channels. A constant injection or 
suction velocity is used. With these assumptions, the Navier-Stokes 
equations for two-dimensional steady state incompressible laminar flow 
are written as 


9u , * 9u 1 9p , 

u-^ + v*^ = - - ^ + u 
9z 9y /) 9z 


9^u , 9^ 
" 9z^ 9y 


i] 


(1) 


“az ay %y2 J 


(2) 


and the continuity equation is 


M 4 . n 

9z ay ■ ^ 


( 3 ) 


The boundary conditions used on the channel surfaces are 


u = 0, and v* = 0 
u = 0, and V* = Vq* 
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at y = 0 (solid wall) 
at y = D (porous wall) 
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Figure 1. Schematic diagram of a semi porous channel. 
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Similarity Transformations 


The governing equations can be transformed into total 

differential equations by the use of a dimensionless length 
coordinate, X, and a dimensionless function, f, which automatically 

satisfies the continuity equation. The new variables are defined as 

X » ^ (6) 

f(x) • ^ 0 ) 

V *z 

where V(z) = U*(0) - 

and U*(0) is the average velocity at z =0. Then, the local velocity can 
be expressed by the new variables as follows 


u = - ^ = Vf 
“ D ay 


( 8 ) 


V* = - M = V. f 

3z 0 


(9) 


Use of equations(l) , (2), and (6) to (9) gives 


_ i ^ = V 

p 3z '' 


-2_ (ff. . (f.)2 - ^ f* 


( 10 ) 
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p ax 




■I 


( 11 ) 
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since the right side of equation (11) is a function of X 


only, differentiating with respect to X yields 


3zax 


On taking the derivative with respect to 
and integrating, equation (10) becomes 


( 12 ) 

X, using equation (12), 


RSq [(f)^ - ff"] + f" = - A (13) 

v.*D -.2 

where Re. = — — and A = - — ^ 

0 1/ i/p Qx 

The transformed boundary conditions are written as 


f = f = 0 

at X = 0 (solid wall) 

(14) 

II 

o 

u 

at X = 1 (porous wall) 

(15) 


Solution of Equation 

The differential equation(13) together with the associated 
boundary conditions constitute a nonlinear boundary value problem with the 
parameter ReQ. Since the governing equations are no longer partial 
differential equations, they can be solved numerically by the Runge- 
Kutta integration method. However, the number of boundary conditions is 
not enough to solve the equation(13) . Therefore, for each 
specified Reo, the value of A and f"(0) are guessed to solve the system 
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of the first order differential equations.. Then, at X = 1 (porous 
wall), the calculated values of f and f are compared with the given 
boundary conditions. When both values are not acceptable, other values 
are used until the tolerance is satisfied. When wall Reynolds numbers for 
suction are greater than 13, separation occurs on the solid wall. 
Similarity solutions are valid for wall Reynolds numbers up to 13. Hence, 
the equation is solved for wall Reynolds numbers ranging from - 30 to 13. 

A comparison of the present results with those found in 
the literature [3,4], which used the perturbation method for small 
Reynolds numbers, and a different numerical scheme, verifies the 
accuracy of the numerical solution. 

The velocity profiles for the suction and injection sections are 
shown in Figures 2 and 3. Asymmetric boundary conditions means that the 
velocity profiles are asymmetric. For injection, Reo<0, the location of 
the maximum velocity shifts from the center of channels to the 
solid wall. Thus, injection increases friction at the solid wall, and 
decreases friction at the porous wall. As wall Reynolds numbers 
increase, the degree of shifting is large. However, the general 
shape is changed little. For suction, Roq > 0, the location of the 
peak velocity shifts toward the porous wall. Hence, suction increases 
friction at the porous wall and decreases friction at the solid wall. 
Unlike injection, separation appears around a wall Reynolds number of 13. 
The velocity profiles change considerably with wall Reynolds numbers. 
In general, friction for the semiporous channel is larger than that for the 
impermeable channel. 
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2. Compressible Vapor Flow Analysis 
Steady, compressible, one-dimensional laminar flow in a heat pipe 
is considered. The heat pipe under study is shown schematically in 
Figure 4. For purposes of formulating the mass, momentum, and energy 
equations in one-dimensional form, the velocity is taken to be the 
average velocity, which is approximated by the velocity distributions based 
on the similarity solution of semi porous channels. Shear stress at 
the interface and the momentum and energy factors are similarly calculated 
and shown in Figure 5. The vapor is evaporated at the liquid-vapor 
interface with mass injection rate, mg , per unit area. It is assumed 
that this vapor flows inward with a normal component of velocity only, 
and joins the axial flowing stream. 

Principle Governing Equations 

The principles of conservation of mass, momentum, and energy in 
the axial direction, and the previously mentioned assumption, yield 



(16) 


dz 




80 


(17) 



( 18 ) 
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Figure 4. Heat pipe vapor region. 



(X)d2i 


12 - 


Figure 5. Variations of momentum and enei 



where D is the thickness of the heat pipe and h is the vapor enthalpy. 
The friction factor, F, of the surfaces is written as 


8r 

F = — ^ 


( 19 ) 


The momentum and energy factors, a and p, respectively, are expressed 


as 


a = 



dy 


( 20 ) 



( 21 ) 


The normal velocity of the vapor at the interface is expressed in 
terms of the heat flux and latent heat of vaporization as follows 


V 


0 


* 



( 22 ) 


The fluid in the vapor passage of the heat pipe is assumed to be 
a mixture of liquid and monotonic vapor. Thus, the quality of the vapor 
is considered. The specific volume v and the enthalpy h are expressed as 


V = Vf + X • (Vg - Vf) 


h = hf + X • hfg 


( 23 ) 

( 24 ) 
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The specific volume of the saturated vapor can be approximated by 



(25) 


Also, the temperature and pressure are related by the Clausius- 
Clapyron equation: 


P ■ R j2 


(26) 


where hfg is the enthalpy of saturated liquid, vf and Vg are the specific 
volume of saturated liquid and vapor, respectively, and X is the vapor 
qual ity. 

Formulations of Differential Equations 

Bankston et. al.[5] showed that variation of the momentum and 
energy factors with axial distance and radial Reynolds numbers is very 
small except near the end of the heat pipe. Also, Figure 5 shows both a 
and p are nearly independent of the radial Reynolds numbers except near 
separation. These results • are taken from the calculations for the 
similarity solution of the semi porous channels. Therefore, it is assumed 
that the derivative of a and p with respect to axial distance is equal 
to zero. In addition, assuming Vf and hfg are constant quantities, and 
combining equations (16) thru (26) yields the axial gradients for the 
density, quality, velocity, pressure, and temperature. 
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Density 

The differentiation of equation (25) with respect to z gives 


dZ ■ " P dZ T dZ 


( 27 ) 


Substitution of equation (26) into the equation above yields 


dv_ v„ 

—2 = -2 
dZ P 


RT 

L F 


- 1 


fg 


dZ 


(28) 


The derivative of specific volume in the z-direction is written as 


S-Jm, 


■|| SI • > 


dv_ dVx 

_2 I 

. dZ dZ 


(29) 


With the prescribed assumption for vf, substitution of equation (28) 
into (29) obtains 


dZ 



RT 

I F 


- 1 


fg 


dZ 


(30) 


Also, the derivative of specific volume and density are related 
as follows 


^ _ i_ ^ 

dZ - %2 dZ 

From equations (30) and (31), the 
differential form can be written as 


(31) 


expression of density in 
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dP 
dZ J 


( 32 ) 


Qua! i ty 

From equation (9), the quality gradient is expressed as 


dx V 

dZ “ (Vg - v^) 


L dZ ^ P 


1 - 


RT 

''fg ^ 


dP 
dZ J 


Substitution of equation (31), equations of the 
mass, and momentum into equation (33) yields 


(33) 


conservation of 


dx 

dZ 


,2.. V 


(Vg - V,)V‘ 




BI_ 11 ^ . FV‘ 
h^g JJ dz “ 8Dva 


(34) 


(Vg - v^) ' DV 

Velocity 

The equation of the conservation of mass can be expressed as follows 


vm. 


dZ ~ 0 V dZ 


(35) 


Substitution of equation(30) into the expression above gives 
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( 36 ) 


^ ^ ~ ^ ^ ^ f BL 1 1 ^ 

dZ ‘ D V dz P V I h^g “ ^ J dZ 

Pressure 

With the previously mentioned assumptions, substitution of the 
energy equation (18) into the momentum equation yields 


^ ao dh °*"o 
dZ " ^ dZ ■ ^VD 



(37) 


In order to obtain an expression for the derivative of the 
vapor enthalpy, differentiation of equation (24) is taken with respect 
to Z, and equation (26) is substituted: 


^ h ^ 4. . I? ^ 

dZ ■ "fg dZ h^g P dZ 


(38) 

enthalpy and 


After the expressions for the derivative of the 
quality are substituted into equation (36), the pressure gradient can be 
expressed as follows 


dZ 


fo f, 

■ /?vo r 


/?v0 i^*’fg'^ 


V!f 


[ho-h*^j,-+- 2 - JJ- J 8 D - ‘ 


“ P „2 


a , ''fq^ , ''q [. RT 1 a '^pfR 

P ” v2 P "fgJ P "fg ^ 


(39) 
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Temperature 

The temperature gradient is derived from the Clausius-Clapyron 
equation as follows 


^ riL . ^ 

dZ ■ h^gP dZ 


(40) 


Description of Numerical Scheme 

The five dependent variables, i.e. density, quality, 
velocity, pressure, and temperature, are coupled by the differential 
equations (32), (34), (36), (39), and (40) which are first order, 
nonlinear, and ordinary differential equations. A computer code has 
been written to solve these equations numerically and simultaneously 
using the Runge-Kutta integral method, which needs only proper 
boundary conditions for the first step of integration. 

Since the velocity term appears at the denominator of equation (39), 
the physical boundary condition for velocity can not be used directly. To 
avoid this problem, proper boundary conditions at the upstream end 
of the evaporator are determined to initiate all calculations as 
follows: a new boundary condition for the velocity is determined at a 
short distance away from the beginning of the evaporator, assuming 
incompressible and saturated vapor flow at the temperature corresponding 
to the heat pipe operating temperature. Also, the saturation pressure 
is used for the pressure. Since the velocity is small at this point 
close to the upstream end of the evaporator, boundary conditions 
determined this way are realistic. The differential equations are solved 
by using these boundary conditions. 
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Results and Discussion 


To investigate various effects, two types of heat pipes with sodium 
as the working fluid were selected for the analysis. One has an 
adiabatic section, the other does not. The dimensions of the heat pipes 
are shown in Table 1. Also, two different operating temperatures are 
chosen such that, for low operating temperatures, the pressure drop along 
the heat pipe is relatively large at the same heat flux. For each 
operating temperature, several cases of heat pipe operation with 
various heat fluxes are considered. Since the similarity solutions 
used are valid only while wall Reynolds numbers are less than 14, the 
maximum heat flux on the condenser is chosen to satisfy this condition. 


Dimension 

773 OK 

8 O 8 OK 

Le (cm) 

8.0 

20.0 

La 

0.0 

10.0 

Lc 

22.0 

30.0 

D 

1.33 

1.35 

W 

1.33 

1.72 

Ao(cm2) 

10.64 

34.4 


Table 1. Primary Information on Sodium Heat Pipes. 
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Figure 6 shows axial variations of vapor temperature, 
pressure, velocity, and density which are obtained from equations 
(32), (34), (36), (39), and (40). The operating temperature of 773 K, 
and the uniformly distributed heat transfer rate of 200 watts are used. 
Momentum equation (17) implies that the variation of pressure in the axial 
direction depends on the relation between the contributions from 
inertia and friction. Thus, the pressure of the vapor in the evaporator 
falls sharply along the vapor passage due to friction and the 
acceleration of the flow due to the injection of mass. The 
corresponding temperature also drops sharply about 5.5 K. Since the 
density of sodium vapor is relatively low at the low operating 
temperature, the velocity is correspondingly large to transfer the 
required mass. The results show the maximum Mach number reaches 0.3 at 
the exit of the evaporator. 

In the condenser, the extraction of mass leads to a deceleration of 
the vapor flow so that pressure increases, while friction reduces the 
pressure. Therefore, the pressure can not recover completely due to 
friction loss. Even though temperature recovery is achieved in 
the condenser, the temperature at the end of the condenser is about 4.2 
K less than that at the beginning of the evaporator. Therefore, the 
vapor passage in the heat pipe is not isothermal. 

Figures 7, 8, and 9 show variations of temperature, pressure, and 
Mach number, respectively, corresponding to five different heat fluxes. 
At the same operating temperature of 773 K, a larger heat input leads to 
greater pressure and temperature drops, and to a higher Mach number. When 
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Figure 6. Axial variations of vapor temperature, pressure, velocity, and density at operating temperature 
of 773 k. 
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Figure 7. Axial variations of vapor temperatures for various heat inputs. 
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the Mach number exceeds about .2, the variatipn of the Mach number is 
large due to expansion of the vapor. For heat inputs of 50 and 100 
watts, the vapor can be assumed to be isothermal. 

Figure 10 shows that the possible maximum heat transfer is limited 
to about 351 watts because the Mach number at the end of the 
evaporator approaches the sonic velocity. 

Figure 11 presents axial variations of vapor temperature, 
pressure, velocity, and density at the operating temperature of 808 K. For 
this case, the operating temperature is 35 K higher, but pressure and 
density are twice as high than for 773 K. Thus, the velocity is small 
and the maximum Mach number is only about 0.2. This leads to a small 
pressure drop, and the corresponding temperature, drop is about 2.8 K. 
Even though the heat input is twice as much as that for 773 K, a 
comparison of temperature drops implies that pressure and temperature 
drops mainly depend on the operating temperature. 

In the adiabatic section, there is no mass injection or 
extraction. Hence, pressure simply decreases due to friction, and the 
velocity increases very little. The contribution of the adiabatic section 
is the addition of pressure and temperature drops such that for the 
high Mach number at the exit of the evaporator, the maximum heat 
transfer rate is reduced considerably. Figure 12 shows this 
phenomenon. The heat input of 865 watts results in a Mach number of 0.7 at 
the exit of the evaporator. Then, the vapor is accelerated in the 
adiabatic section, so the velocity of vapor approaches the sonic velocity 
at the end of this section. 
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Figure 11. Axial variations of vapor temperature, pressure, velocity, and density at operating 
temperature of 808 k. 
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In the condenser section, as expected, , the pressure and 
temperature drops are recovered, but the degree of recovery is small due 
to small vapor velocity and large friction loss in the longer condenser. 
Numerical results show that the temperature recovery is only 0.6 K. 

Variations of temperature, pressure, and Mach number corresponding 
to six different heat fluxes at the operating temperature of 808 K 
are presented in Figures 13, 14, and 15, respectively. For heat input up 
to 200 watts, pressure recovery is not observed. This implies that the 
pressure drop due to friction in the condenser may exceed the value of the 
inertial contribution. However, the total pressure drops are small, so the 
vapor can be assumed to be isothermal. 

Figure 16 shows the distributions of the nonuniform heat flux on 
the evaporator and condenser. Axial variations of vapor temperature, 
pressure, velocity, and density corresponding to the distribution in Figure 
16. a are presented in Figure 17. At a short distance away from the 
beginning of the condenser, the slight heat extraction suggests that the 
axial variation of the vapor velocity is small. Thus, the friction 
effect in this region is dominant. The minimum temperature and pressure 
appear in the condenser instead of at the exit of the evaporator. As 
the heat extractions increase, the absolute value of axial variation of 
velocity increases and the vapor velocity decreases. The pressure 
drops are recovered gradually. In Figure 18, corresponding to the 
distribution in Figure 16. b, the temperature and pressure drops are 
recovered immediately from the beginning of the condenser due to the 
large heat extraction. As expected, axial variations of 
temperature, pressure, velocity, and density depend on the distributions 
of heat input. 
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Figure 13. Axial variations of vapor temperatures for various heat inputs. 
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Figure 16. Distribution of heat input on the evaporator and condenser. 
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Figure 17. Axial variations of vapor temperature, pressure, velocity, and density for nonuniform 
heat input in figure 16a. 
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Figure 18. Axial variations of vapor temperature, pressure, velocity, and density for nonuniform 
heat input in figure 16b. 



In summary of the results, the degree of pressure recovery 
mainly depends on the heat flux, operating temperature, and length 
of the condenser. The higher the heat flux, the lower the operating 
temperature, and the shorter the condenser, the greater the pressure 
recovery. However, for this case, the pressure drop in the evaporator can 
be large enough that the entire heat pipe is not isothermal. 

Since experimental data is not available, quantitative comparison 
can not be achieved. However, comparison of the general behavior of the 
present results with published data[2,6,7] gives qualitative agreement. 
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LIST OF SYMBOLS 


Aq 

Cpf 

D 

f 

F 

h 

hf 

hfg 

ho 

*-e 

*-a 

Lc 

ffio 

P 

Q 

R 

Reo 

T 

u 

U*(0) 

V 

Vf 


area of porous wall of injection or suction section 

constant pressure specific heat of liquid 

thickness of vapor passage 

dimensionless stream function 

friction factor 

enthalpy of vapor 

enthalpy of saturated liquid 

latent heat of vaporization 

enthalpy of vapor at interface 

length of evaporator 

length of adiabatic section 

length of condenser 

mass injection or suction rate per unit area 

pressure 

quality of vapor 

gas constant 

wall Reynolds number 

temperature 

z-component of velocity 
average velocity at z = 0 
specific volume of vapor 
specific volume of saturated liquid 
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Vg specific volume of saturated vapor 

v*o velocity at porous wall 

V* y-component of velocity 

V average velocity at cross section 

W width of vapor passage 

X vapor quality 

y coordinate direction normal to axial direction 

Z coordinate in axial direction 

Greek Symbol 

OL momentum factor 

p energy factor 

X dimensionless length coordinate 

V kinetic viscosity 

p density of vapor 

Pq density of vapor at interface 

Ty shear stress 

stream function 
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